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Abstract 

We present a hybrid Brownian dynamics / Monte Carlo algorithm for simulating solutions of 
highly entangled semiflexible polymers or filaments. The algorithm combines a Brownian dynamics 
time-stepping approach with an efficient scheme for rejecting moves that cause chains to cross or 
that lead to excluded volume overlaps. The algorithm allows simulation of the limit of infinitely 
thin but uncrossable threads, and is suitable for simulating the conditions obtained in experiments 
on solutions of long actin protein filaments. 
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I. INTRODUCTION 



Experimental and theoretical studies of solutions of entangled semiflexible polymers have 
probed concentration regimes that are difficult to simulate by standard methods. In solutions 
of actin protein filaments, the typical filament length L and persistence length L p are both 
of order 10/xm, but the chain diameter d ~ 8 nm is about 10 3 times smaller. For such 
solutions at concentrations near the isotropic-nematic transition, the actin volume fraction 
is of order 10~ 3 , cL 3 ~ 10 3 — 10 4 (where c is a number of chains per unit volume), and 
the typical distance between chains is of order 0.1/im. These conditions may be fruitfully 
idealized in theoretical work by considering the limit d = of uncrossable but infinitely thin 
threads. A simulation of these conditions with a conventional molecular dynamics (MD) 
or Brownian dynamics (BD) model of polymers as necklaces of nearly-tangent repulsive 
spherical beads jl] would require a large number N ~ L/d ~ 10 3 beads per chain, and an 
extremely small time step in order to resolve both the strong short range repulsion and the 
bending forces needed to maintain angular correlations over 1000 beads. Here, we present 
a hybrid Brownian Dynamics (BD) / Monte Carlo (MC) algorithm for simulation of such 
systems of very thin threads, in which chains are represented as sequences of a much smaller 
number of uncrossable rods. The algorithm is designed to allow efficient simulation of the 
idealized limit of infinitely thin but uncrossable threads, as well as of chains with a small 
but nonzero steric diameter. 

The article is organized as follows. Sec. |H] is a description of our algorithm for simu- 
lating solutions of uncrossable bead-rod chains. Sec. IIHI is a description of the geometrical 
algorithms used to efficiently detect chain crossings and to calculate the distance of closest 
approach between rods. Sec. US presents results of tests of the validity and computational 
cost of the algorithm. Sec. is a discussion of the theory underlying the algorithm. Con- 
clusions are summarized in Sec. IVII 



II. SIMULATION ALGORITHM 



Our simulation method combines a Brownian dynamics time-stepping algorithm with a 
scheme for rejecting moves that cause chains to cross or overlap. At each step, a Brownian 
dynamics (BD) algorithm that has been used previously to describe non-interacting wormlike 
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chains is used to generate a trial move for a randomly chosen chain. To simulate a solution 
of infinitely thin but uncrossable threads, we reject all trial moves that cause one chain to 
cross through another, and accept all others. 



A. Single Chain Brownian Dynamics 

The simulations use a discretized model of wormlike chains in which each chain in a 
solution is represented as a set of N + 1 beads, which act as point sources of friction, 
connected by N inextensible rods, each of which is constrained to have a fixed length a. Let 
R M be the position of bead fi of a chain, with \x = 1, . . . , N + 1, and qj = R»+i — Rj be 
the vector associated with rod i, which is constrained to have a fixed length |q$| = a, for 
all i = 1, . . . , N. Let Uj = qj/|q*| be a unit tangent vector for rod i. The bending energy, 
denoted by Uq, is 

N 

u = — VV-u^i , (i) 

a ' 

i=2 

where k is the bending rigidity of the chain. The persistence length is L p = n/ksT. This 
model approximates the bending energy of a continuous wormlike chain in the limit a <C L p 
and N > 1. 

To generate trial moves for individual chains, we use a Brownian dynamics algorithm for 
chains with constrained rod lengths that has used previously to study dilute solutions of 
both semiflexible and flexible chains 8]. In this algorithm, changes in bead 

positions are calculated from a bead velocity V M of the form 



9U Q ^ ^ , 



(2) 



Here, H M is a mobility tensor for bead /i, % is a tension in rod i, and r)^ is a random 
Langevin force. F™ et is a metric correction force that is required in order for this algorithm 



to yield the correct equilibrium distribution 



We use an anisotropic bead mobility tensor of the form 



.71 



H M = -^u^u^ + -^-(1 - u^u^) , (3) 

Sll S-L 



where u M is an approximation to the unit tangent vector at bead i, and Oi and C4_ are friction 



coefficients per unit length for parallel and perpendicular motion, respectively. p| . The unit 
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tangent u M is approximated by a centered difference u M = (R M +i — R M _i)/|R M +i — R^-il for 
beads \i = 2, . . . , N, and by ui = ui and Ujv+i = for the end beads. 

The random force rj^ is generated by the procedure described by Morse and Montesi 
et al. 0] for generating "geometrically projected" random forces. At the beginning of each 
time step, an unprojected random force 77^ for each bead is generated from a distribution 
with a vanishing mean, (77^) = and a variance (rf^rjl) = 2/«bTH~ 1 <5 a1 j,/ At. The force r?^ in 
Eq. (j2J) is a projected force of the form 77^ = 77^ — r^u^ — i^-iu^x in which the quantities 
fji, . . . , ?7jv are calculated by requiring that = Uj • {r] i+1 — T7J for alH = 1, . . . , N. 

The tensions are calculated, after calculation of the bending, metric, and 

random forces, by requiring that the rod lengths all maintain constant length, or that = 
Uj • (Vj + i — Vj) for all % = 1, . . . , TV. This yields a tridiagonal set of N equations, which 
must be solved every time step. 



Each time step for a sing 
inally proposed by Fixman 



e chain is generated by a mid-step integration algorithm orig- 
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lOj]. At the beginning of each time step, projected random 



forces are generated. Predicted mid-step bead positions are then calculated as 

R (l/2) = R (0) + V (0) At/2 (4) 

where R^ represents an initial bead position and represent a bead velocity calculated 
using the initial bead positions. Final bead positions are calculated as 

R« = Rf +V^) A t (5) 

where V^ 1 1//2 ' ) is a bead velocity computed using bending and metric forces, mobilities, and 
tensions that are re-calculated using the mid-step bead positions, but using the same random 
force 77^ as that used to calculate the mid-step positions. 



B. Detecting chain crossing 

To simulate a solution of infinitely thin but uncrossable threads, we combine this Brow- 
nian dynamics algorithm with a scheme for rejecting moves that cause chains to cross. At 
every timestep, a trial move is generated for a randomly chosen polymer chain using the 
Brownian dynamics algorithm discussed above. The trial move is rejected, and the chain 
position is unchanged, if the trial move would cause any of the N rods of the chosen chain to 
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cut through any rod of any other chain. Whether or not the move is accepted, another chain 
is then chosen at random, and the process is repeated. To accurately simulate Brownian 
motion, the time step At used in the Brownian dynamics algorithm must be chosen to be 
small enough so that almost all of the moves are accepted. In a system containing M chains, 
a time At is taken to have elapsed every M attempted chain moves. 

The usefulness of the algorithm relies upon the availability of an efficient method for 
detecting when a trial move causes one chain to cut through another. Consider an attempted 
move of a chain a. To detect whether this attempted move would cause chain a to cut 
through any other chain, we consider each of the N rods of chain a to sweep out a surface 
over the course of a time step. We reject the move of the chain if the surface swept out 
by any rod of a is intersected by any rod of any other polymer /3, as show schematically in 
Figure HJ 

Consider the test for whether a move of rod i of the moving chain a will cause it to cut 
through a rod j of a stationary chain (3. The line segment corresponding to the stationary 
rod (rod j of (3) may be described by a parametric equation 

r'O') = c' + s'q'. (6) 

where c' is the center of the rod, q' is the vector connecting its ends, which is of length 
|q'| = a, and s' is a contour length parameter with a domain —1/2 < s' < 1/2. The surface 
swept out by the moving rod (rod % of a) during a trial move can be represented by a two 
parameter equation, 

r{s,t) =c{t) + sq{t), (7) 

where — 1/2 < s < 1/2, t is a dimensionless time variable that ranges from to 1 during 
the time step, and c(t) and q(t) are the center position and end-to-end vector of the rod as 
functions of t. Thus, r(s, 0) and r(s, 1) represent rod conformations at the beginning and 
end of the timestep, respectively. To define the surface swept out at intermediate times, we 
take c(t) and q(t) to be linear functions of t, 

r(t) = c + t Ac 

q(t) = q + t Aq (8) 

where c = c(0), qo = q(0), Ac = c(l) — c(0), and Aq = q(l) — q(0). This is equivalent 
to assuming that the beads at both ends of the rod move with constant velocity from initial 
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FIG. 1: Intersection of a stationary rod with the surface swept out by a moving rod 



to final position, and that the rod is always a straight line between these two beads. The 
length of the rod is not exactly conserved at intermediate times, but the surfaces swept out 
by neighboring rods within a chain meet without gaps or overlaps along the trajectory of 
the shared bead. 

For the surface r(s,t) swept out by the moving rod to be intersected by a stationary 
rod r'(s'), the line segments r'(s') and r(s, t) must intersect at some time t = tj, such that 
< ti < 1, at a point of intersection for which —1/2 < s < 1/2 and —1/2 < s' < 1/2. An 
efficient algorithm for determining whether such an intersection occurs for any specified pair 
of rods is described in Sec. IIHI 
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C. Rod- Rod Interactions 



The algorithm for infinitely thin threads can be generalized to allow for intermolecular 
interactions. We consider a model in which the total potential energy is a sum of a bending 
energy Uq and an interaction potential Uint- The interaction energy is taken to be a pairwise 
additive sum of interactions between rods, 

U iat = Y f U{\d ij \), (9) 

i>j 

where the sum is taken over all pairs of rods, and where the interaction energy Z7(|dj_j|) 
between rods i and j is a function of the distance of closest approach between those 
two rods. An algorithm for calculating distance of closest approach between two finite line 
segments, and a discussion of how this distance is defined, is given in Sec. IHII 

In an algorithm in which trial moves are generated using a Brownian dynamics algorithm 
for non- interacting chains with potential energy Uq, trial moves that do not cause chains to 
cross are accepted or rejected according to a Metropolis acceptance criterion that is based 
on the change in interaction energy: A trial move is accepted if it causes Ui nt to decrease, 
and is accepted with a probability 

P acc = exp[-AU mt /{k B T)} (10) 

if it leads to an increase AUi n t > in Ui n t. In the case of a hard-core interaction, a move 
that does not cause chain crossing is rejected if it causes the distance |dy| between any two 
rods to fall below a specified hard-core diameter, and accepted otherwise. Note that only the 
interaction energy is used in this acceptance rule, and not the change in the intramolecular 
bending potential, because the effect of the bending potential is already fully accounted for 
by the Brownian dynamics algorithm used to generate trial moves. The rationale underlying 
the use of such a Metropolis sampling is discussed in Sec. El 

Each step of the simulation thus involves the following sequence of operations. A trial 
move is generated for a randomly chosen chain. A check is performed for every rod on 
that chain to see if the trial move would cause intersection of that rod with any of the 
neighboring rods, which are identified by use of a Verlet neighbor list. The move is rejected 
if any intersections are found. If no intersections are found, and Ui n t ^ 0, the change in 
the interaction energy of each rod with all nearby neighbors is calculated, and the move is 
accepted or rejected based upon the Metropolis criterion described above. 
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111 ]. Verlet neighbor lists are used to reduce 



As in simulations of systems of point particles 
the number of rods that must be checked for intersections and for which interaction energy 
must be calculated. The neighbor list for each rod in the system contains all other rods for 
which the distance of closest approach was less than a Verlet radius when the neighbor 
list was constructed. The neighbor list is reconstructed whenever any bead in the simulation 
is found to have moved a distance greater than since the last update of the list. The 
neighbor list is constructed with the use of a linked cell list of rods, by a procedure closely 
analogous to that used for point particles. Optimal values of the Verlet radius for each set 
of parameters have been determined by extensive trial and error. 



D. Preparation of initial states 

Many of the simulations for which we have used our algorithm focus on the dynamics and 
stress relaxation in entangled solutions over relatively short time scales. These simulations 
are often run for times shorter than a reptation time due to computational limitations on our 
ability to simulate very crowded solutions for long times. For such short simulations to yield 
valid results, they must start from initial configurations that are already representative of a 
thermal equilibrium ensemble for the solution. Such simulations thus rely upon our ability to 
efficiently generate equilibrated initial states for systems that would be difficult or impossible 
to equilibrate by running the algorithm described above, which we use to study dynamical 
properties. 

In the case of infinitely thin but uncrossable chains, with Ui n t = 0, the prohibition 
on chain crossing has an enormous effect upon dynamics of the system, but has no effect 
upon the equilibrium probability distribution, since it changes neither the set of allowed 
microstates nor the energy of any microstate. The equilibrium distribution for a solution of 
such uncrossable chains is thus identical to that of a solution containing an equal number 
of non-interacting "phantom" chains. This distribution is characterized by random chain 
positions and orientations, and by a Boltzmann distribution 

e (L p /a)cos 8ij 

P(cos 6ij) = == — — , r , , — — , (11) 

Jo rf%sin^e(^/»)co S 0y' v > 

for the cosine cos 9ij = u; • Uj of the angle between neighboring rods j — % ± 1 within a chain. 
We may thus generate initial configurations for such a solution by a method used previously 
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to generate initial configurations for dilute solutions, in which chains are grown by a process 
that yields an equilibrium distribution of chain conformations. The first bead of each chain is 
placed at random within the simulation cell, and the first rod is given a random orientation. 
The orientation of each subsequent rod relative to the previous one is chosen such that the 
cosine of the angle between these two rods is chosen from the Boltzmann distribution given 
in equation (fTTj) . 

When Ui n t 7^ 0, inter-molecular interactions lead to non-trivial correlations. In this 
case, we start with an initial configuration for a non-interacting system, generated by the 
method described above, and then equilibrate the system by a slithering- snake Monte-Carlo 
algorithm. In each attempted move of this algorithm, a rod is removed from one end (the 
"tail") of a randomly chosen polymer and reattached to the other end (the "head"). The 
orientation of the new rod relative to that of the pre-existing end rod is chosen from the 
equilibrium distribution given in where cosO is the angle of the last joint in the chain. 

The resulting attempted Monte-Carlo move is accepted or rejected based on a Metropolis 
criterion based on the resulting change in interaction energy Ui nt , as given in Eq. ()10|) for 
the acceptance probability of a move that changes the interactions by an amount Ui nt > 0. 
Monte Carlo equilibration is continued until every chain in the simulation box has reptated 
many times its own length. In the concentration range of interest, where the system is 
highly entangled but semi-dilute, and in which the equilibrium phase is isotropic rather 
than a nematic liquid crystal, almost all such Monte Carlo moves are accepted, and the 
system can be equilibrated quite rapidly. 

III. GEOMETRICAL ALGORITHMS 
A. Intersection 

To determine whether two line segments intersect during a move, we first determine 
whether the two infinite lines containing these segments intersect during the time < t < 1 
of interest. If an intersection time < tj < 1 is found, we then calculate the contour 
variables s and s' at the point of intersection, and determine if these lie within the two line 
segments that define the rod. 

The efficiency of our algorithm relies upon the existence of efficient algorithms for de- 
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tecting intersections between rods and calculating distances of closest approach between line 
segments. 

An intersection between a moving line r(s, t) and a stationary non-parallel line r'(s') that 
pass through rod centers c(t) and c', respectively, can occur at a time tj if the vector c(tj)—c 
has a vanishing projection onto a vector q(ir) x q' that is perpendicular to both lines, i.e., 
if 

0=[c(t / )-c']-[q(t/)xq / ] (12) 
Substituting Eq. (JHJ) for c(i) and q(i) into this condition yields a quadratic equation 

at] + + c = . (13) 

for the intersection time tj of the two infinite lines, in which 

a = Ac • (Aq x q') 

b = Ac • (q x q') + [co - c'] • (Aq x q') 

c = [c -c']-(q xq') (14) 

If a real solution is not found in the domain < tj < 1, then these two lines do not intersect 
during the time step of interest. 

If the infinite lines containing the rods do intersect at a time < ti < 1, we must then 
calculate the coordinates s and s' at the point of intersection, where r(s, tj) = r(s'). The 
two rods intersect only if these coordinates satisfy —1/2 < s < 1/2 and —1/2 < s' < 1/2. 
We calculate the coordinates s and s' as a special case of the algorithm for determining the 
point of closest approach of two lines, which is discussed below. 

Ideally, these geometrical checks should be performed for every rod of the moving chain 
a and every other rod in the simulation box, including other rods on chain a. The algorithm 
described above, however, is designed to check only for intersections of rods on a moving 
chain with rods on a different stationary chain. It is possible to generalize the algorithm 
so as to identify intersections between two moving rods. The generalization requires that 
the quadratic equation for tj, equation (JHJ), be replaced by a cubic equation. However, 
intersections between different parts of the same chain are expected to be rare in the systems 
of rodlike chains upon which we have focused and so, for simplicity, we have not implemented 
a check for such self-intersections. 
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B. Distance of Closest Approach 



Consider two rods that correspond to line segments — 1/2 < s < 1/2 and — 1/2 < s < 1/2 
of the lines 

r(s) = c + sq 
r'(s') = c + s'q! . (15) 

The distance of closest approach (DCA) between these two rods is the minimum magnitude 
\d(s, s')\ of the separation 

d(s, s') = r(s) - r(s') (16) 

for — 1/2 < s < 1/2 and —1/2 < s' < 1/2. This minimum may be obtained when s and/or 
s' are equal to ±1/2, in which case the distance of closest approach is the distance between 
the end of one rod and a point along the length of the other, or the distance between two 
rod ends. 

The first step in the calculation of this minimum distance is the identification of the 
coordinates s and s' at the point of closest approach of the infinite lines that contain the 
two rods. To identify this point of closest approach, we require that 

= d\d\ 2 /ds = 2q-d(s,s') 

= d\d\ 2 /ds' = 2q'-d(s,s') (17) 

Note that this is equivalent to the geometrical requirement that the vector d(s, s') at the 
point of closest approach be perpendicular to both q and q'. Substituting Eqs. (fT5j) and 
(Unj) into these conditions yields a pair of linear equations 

as-bs' = d (18) 
-bs + cs' = e (19) 

in which a = |q| 2 , b = q • q', c = |q'| 2 , d 
solution 

s - 
s' -- 



q • [c' — c] , and e = q' ■ [c — c'] , which have the 



be + cd 
ac — b 2 

^ . (20) 
ac — b 2 
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As part of the check for the intersection of two rods, we use Eq. (|2U|) to calculate the 
coordinates s and s' for pairs of lines r(s) = r(s, ti) and r'(s') that have already been shown 
to intersect at a time ti, for < tj < 1. In this case the point of closest approach is actually 
the point of intersection. 

The distance of closest approach between the non-parallel rods is the same as that for 
the corresponding infinite lines if and only if the coordinates s and s' at the point of closest 
approach of the lines satisfy the inequalities —1/2 < s,s' < 1/2. Otherwise, the closest 
approach between rods may correspond to the separation between an end of one rod and a 
point along the length of the other (e.g., s = 1/2 and —1/2 < s' < 1/2), or between two 
rod ends (e.g., s = 1/2 and s' = —1/2). To compute the shortest distance between rods, 
we first use Eq. (}20|) to calculate the coordinates s and s' of the point of closest approach 
for the infinite lines. If both s and s' lie between —1/2 and 1/2, we accept these values and 
calculate the distance \d(s, s')\. If not, we must consider one of two possible types of special 
cases. 

The first special case occurs when one of the two coordinates s and s' at the point of 
closest approach of two infinite lines lies within the range [-1/2,1/2], and the other does not. 
As an example, consider the case s > 1/2 and — 1/2 < s' < 1/2. In this example, the distance 
of closest approach will be between the rod end s — 1/2 and some point —1/2 < s' < 1/2. In 
this example, we may identify the minimum |d(l/2, s')| for unrestricted s' by solving linear 
equation (JT§j) for s' while taking s = 1/2. If the resulting value of s' lies outside [—1/2, 1/2], 
then the DCA between rods is obtained by setting s' to the value s — ±1/2 of the nearest 
end, giving a DCA equal to the distance |d(l/2, ±1/2) | between two rod ends. 

The other type of special case occurs when neither of the coordinates at the point of 
closest approach of the infinite lines lies within [-1/2,1/2]. As an example, we consider 
s > 1/2 and s' > 1/2. In this case we must locate the minima of both |d(l/2,s')| 2 for 
— 1/2 < s < 1/2 and of |d(s, 1 / 2) | 2 for —1/2 < s < 1/2 by the procedure outlined above, 
either or both of which may yield s = s' = 1/2, and accept whichever yields a smaller 
distance. 

The above reasoning applies only to non-parallel rods. Two rods are parallel if and only 
if the denominator in Eq. QUI vanishes, i.e., if ac — b 2 = 0. For parallel rods, if \d\ < \a + b\/2, 
then the DCA is the magnitude of the projection of c — c' onto the plane perpendicular to 
q, which is given by |c' — c — qd/a\. If \d\ > \a + b\/2, then the DCA is the distance between 
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FIG. 2: Probability distribution for cosine of the angle between neighboring rods of the same chain. 
cL 3 = 4000, N = 40, L p = L, G = (±. Legends indicate At in units of r® od . Values in parenthesis 
represent percentage of trial moves accepted. 

two rod ends. 

IV. VALIDATION AND RESULTS 
A. Equilibrium Properties 

To show that our algorithm simulates equilibrated solutions of infinitely thin chains, we 
begin by analyzing results for two equilibrium properties for which analytical predictions 
are available. 

Figure 121 shows simulation results for a histogram of the distribution of angles between 
neighbouring rods. The angular distributions obtained for dilute and concentrated solutions 
agree with each other and with equation (jllj) within our statistical error. 

Figure El plots a radial distribution function obtained from simulations of entangled semi- 
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FIG. 3: Radial distribution function for rods used to discretize polymer chains. cL 3 = 4000, 
N = 40, L p = L, 0| = Legends indicate At in units of r r ° od . Values in parentheses represent 
percentage of trial moves accepted. 

flexible polymers. The radial distribution function P(r) for a test rod i is defined here such 
that P(r)dr is the probability per unit length of the test rod of finding another rod j for 
which the points of closest approach of the lines containing the two rods lies within the 
line segments represented by the rods, and for which the distance of closest approach lies 
between r and r + dr. We show in appendix El that P(r) for a solution of infinitely thin, 
randomly distributed chains is given by 

P(r) = \p (21) 

The radial distribution functions obtained from our simulations agree with that predicted 
by equation (|21j) to within statistical error. This agreement extends to very small values 
of r, comparable to the distance moved by the rod within a single timestep At, for which 
we might have expected our rejection scheme for preventing chain crossing to distort the 
equilibrium distribution. 
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B. Convergence studies 



Results presented in Sec. IIVAI clearly show that statistical equilibrium properties are well 
converged with respect to the integration timestep At. However, this does not necessarily 
imply that dynamic properties are simulated with equal accuracy. 

We expect dynamical properties to be sensitive to the fraction of trial moves that are 
rejected, which increases with increasing time step At. The dependence of the fraction 
of moves rejected on At and other parameters can be understood by the following scaling 
argument: The probability P re j ec t that a chain of N rods and length L = Na will intersect 
another chain during a trial move is roughly P re j ec t ~ P(r)NaAr±, where P(r) ~ p is the 
radial distribution function defined above (the number of other rods with a specified distance 
of closest approach per unit length of the test chain), Ar_|_ is a typical magnitude for the 
transverse displacement of any rod of the chain, and aAr± is comparable to the average 
area swept out by one rod. The r.m.s. transverse displacement of an individual bead or rod 
within a single time step of a discretized chain is of order Ar± ~ ^/D^At of a free bead, 
where = kT/((^a) is the transverse diffusivity of a single bead. By combining these 
expressions, we predict that the percentage of rejected moves scales as 

P reject ~ AcL 3 ^NAt/r? od (22) 

where r® od = (±L 3 /(72kT) is the rotational diffusion time of a rigid rod of length L in 
dilute solution, and A is a universal numerical prefactor. This relation is confirmed by the 
data shown in Figure EJ where Preject is shown to be a linear function of cL 3 a/ N At / r® od in 
systems with several values of cL 3 and N. The dotted line in this figure yields a prefactor 
A ~ 0.078. 

Figure |H1 shows the effect of varying At on a quantity (A<i 2 (t)) that is approximately the 
mean square displacement of the middle bead of a chain along directions transverse to the 
tube. The quantity (A<i 2 (t)) is defined to be the square of the distance of closest approach 
between the position of the the middle bead of a given polymer chain at a time r + t and the 
contour of the same chain at an earlier time r, as shown in Figure 03 In a tightly entangled 
solution, the magnitude of the plateau in this quantity is a measure of the width of the 
tube to which the polymer is confined, and will be analyzed in detail elsewhere. Clearly, the 
data is well converged for all t over the range of timestep values shown here, for which the 
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FIG. 5: Schematic diagram defining (Ad 2 (t)) for a semiflexible beadrod chain with N = 6 

rejection ratio is of order 10 %. Similar results were observed for other dynamic properties. 

Examination of this and other tests of the convergence with respect to time step led us 
to conclude that diffusion in tightly entangled solutions could be accurately represented, to 
within the statistical errors of our simulations, by choosing a time step so as to accept at 
least 90 % of all trial moves. In highly entangled solutions, with cL 3 > 250 our choice of 
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FIG. 6: Convergence of (Ad 2 (t)} with respect to integration timestep. Legends indicate values of 
Ai/r° orf . Values in parenthesis represent the percentage of trial moves accepted. L p = L, G = 
N = 20 for all data shown above. 



time step is controlled by this criterion, rather than by the need to resolve the fluctuations 
of individual chains. Combining this criterion with Eq. yields a simple prescription for 
calculating At. 

Figure [7| shows the effect of changes in the number N of rods per chain on (Ad 2 (t)) in 
solutions with cL 3 = 1000. The behavior of (Ad 2 (t)) at very early times is, of course, very 
sensitive to chain discretization. At longer times, corresponding to the plateau in (Ad 2 (t)), 
the data is less strongly affected by changes in N. 

C. Computational requirements 

Table I compares simulation times and memory requirements from four different simula- 
tions of solutions of semiflexible chains: (1) A simulation of isolated semiflexible chains using 
an algorithm developed previously for simulating dilute solutions, in which only the confor- 
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FIG. 7: Convergence of (Ad 2 (t)} with respect to chain discretization N. Legends indicate chain 
discretization N. Values in parenthesis represent Ai/r^ orf . At/r® od has been chosen such that 
atleast 90% of trial moves are accepted. L p = L and (ji = £j_ for all data shown above. 

mation of a single chain is stored in memory, which uses the same time stepping algorithm 
subroutines as those used here. (2) A simulation of 6912 phantom chains using our code 
for simulating entangled solutions, but without checking for intersections or interactions. 
The primary difference between this and simulation 1 is that it requires the positions of 
all of the chains to be stored simultaneously in memory. (3) A simulation of an entangled 
solution of 6912 uncrossable but infinitely thin semiflexible threads, in which moves that 
cause chains to cross are rejected. (4) A simulation of corresponding solution of uncrossable 
chains with repulsive excluded volume interaction, with a diameter d/L = 10~ 3 , in which 
we must also calculate the change in interaction energy induced by moves that do not cause 
chains to cross, in order to calculate the Metropolis acceptance criterion. All simulations 
were performed on AMD Athlon MP processors with a CPU speed of 1.67 GHz and a cache 
size of 256 Kb, using a code that was implemented in Fortran 90. 
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Simulation 


Vcrlct 

n/L 


Time 
us /rod step 


Memory 
(MByte) 


1) Single 




0.64 


0.0054 


2) Phantom 




5.12 


81.7 


3) Entangled (d = 0) 


0.0375 


7.61 


96.9 


4) Entangled (d = 10" 3 L) 


0.0375 


8.42 


96.9 



TABLE I: Comparison of memory and simulation times required for four different simulations 
of solution of chains with L p = L and N = 40, with a time step (At/r® d ) = 1.40625 x 10 -9 . 
Simulations (2-4) were of a system of 6912 molecules in a periodic cell of volume chosen to give 
cL 3 = 4000. 

From Table I, we see that the simulation of a solution of phantom chains, in which all of 
the chain conformations are stored but the chains do not interact, is about 8 times slower 
than the equivalent simulation of individual chains. This is presumably because the entire 
problem no longer fits in the processors' cache. Once the cost of this change in memory 
structure is paid, the addition of topological checks slows down the simulation only by a 
factor of about 1.5. The addition of a small steric diameter causes an even smaller further 
increase. The times per rod per time step are quite comparable to the times of a few /xs per 
particle per time step required in MD simulations of bead spring polymers with repulsive 
Lennard Jones interactions at liquid-like densities on the same processors. 

V. THEORETICAL BACKGROUND 

The algorithm presented in Sec. |H] can be viewed either as a Brownian dynamics algo- 
rithm or as a form of Monte Carlo algorithm. The Brownian dynamics viewpoint is necessary 
to justify the use of the algorithm to predict dynamical properties. The Monte Carlo per- 
spective is needed to explain the accuracy attained in results for equilibrium properties, and 
to justify our use of a Metropolis sampling scheme to treat the intermolecular potential. 

To highlight some of the conceptual issues that arise in the discussion of our simulations, 
we have found it useful to also consider stochastic algorithms for a related but simpler 
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problem of one- dimensional (ID) diffusion of a random walker with a coordinate x(t) > 
near a hard wall at x = 0. The requirement that the walker cannot go through the wall is 
analogous to the constraint prohibiting two chains from cutting through each other. This 
constraint may be imposed in the ID diffusion problem by simply rejecting moves that take 
the walker through the wall. 



A. Brownian dynamics 



When designing our algorithm for simulating infinitely thin uncrossable threads, we ini- 
tially thought of it as a Brownian dynamics algorithm, in which the rejection of moves that 
lead to chain crossings was viewed as being analogous to the rejection scheme proposed 
above for simulating one-dimensional diffusion near a reflecting boundary. 

To clarify this point of view, and its limitations, we first recall the theory underlying 
a conventional Brownian dynamics algorithm. In this analysis, one considers a Markov 
jump process for some vector of coordinates X = [Xi, . . . ,X n }. in which the conditional 
probability distribution for any random displacement X — > X + AX depends upon an 
adjustable time step At. The relevant set of coordinates X is the set of all bead positions 
in a simulation of a polymer solution or the coordinate x of the random walker in the ID 
diffusion problem. The probability distribution P(X,t) generated by this discrete process 
may be shown to converge in the limit At — ► to the solution of a corresponding Fokker- 
Planck equation 

dP = d(V t P) d(D l3 P) 

dt dX t OXidXj 1 ' 

if the first and second moments of the random displacement AJYj from a state X obey the 
conditions 

VA X) = *„ m 

n ' At^O At 

Dil( X) = , im (24) 
n ; At^o 2 At v ; 

where Vi(X) is a drift velocity vector and D^^X) is a diffusivity tensor. The BD algorithm 

for non-interacting polymers used in our simulations was designed so as to satisfy these 
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conditions for the drift velocity and diffusivity appropriate to a wormlike chain. 

In our simulations of a system of M polymers, the relevant underlying Markov step is the 
motion of a single randomly chosen chain. The actual time elapsed during the motion of a 
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single chain must be understood to be smaller by a factor of M than the timestep of At used 
in the underlying single-chain algorithm. This trivial rescaling of At is needed to compensate 
for a corresponding reduction of the values of the moments (AXi) and (AXi AX j) in Eq. 
(j21| for coordinates X, and Xj, which correspond to Cartesian components of the positions 
of beads on the same chain: Both of these moments are reduced by a factor of 1/M relative 
to the values obtained in a single-chain simulation because there is only a probability 1/M 
that any particular chain will move during a time step. 

The derivation of the analysis outlined above relies upon the assumption that both Vi(X) 
and Dij(X) are smooth functions of X. This condition is not satisfied by diffusion near a 
reflecting boundary, in which the underlying diffusion equation must either be taken to 
contain values of Vi(X) and Dij(X) that jump discontinuously to zero at a boundary, or 
must be explicitly supplemented in the Fokker-Planck equation by a reflecting boundary 
condition. Several authors have considered how to correctly implement a reflecting boundary 



condition in a stochastic simulation, anc 
from the presence of such a boundary 



proposed methods for minimizing the error arising 



ziposcd i 



The correct probability distribution 



function P(X,t) can actually be obtained in the limit At —>■ from a variety algorithms 
in which the random walker is prevented from penetrating the wall. For example, for a ID 
walker confined by a wall to i > 0, the same probability distribution is obtained in the 
limit At — > either by rejecting moves that yield a trial position x < 0, or by replacing 
trial values of x < by — x. The choice of an algorithm for treating jumps near a reflecting 
boundary does, however, affect the magnitude of the systematic time-discretization error 
produced by a stochastic simulation for At ^ 0. In the simple ID diffusion problem, in 



which each particle jumps a distance of order Ax ~ yDAt per time step, most "plausible" 
rules for treating jumps near a reflecting boundary introduce a large, 0(1) error in P(x) 
over a region of size O (Ax) of the boundary, and produce smaller 0(-\/~Ab) errors far from 
the wall. 

In light of this understanding, we thus initially expected our algorithm to yield a radial 
distribution function P(r) for the distance of closest approach between pairs of uncrossable 
rods that deviates significantly from theoretical predictions for inter-rod distances r of order 
the typical displacement Ax of a single bead over one time step, which is given by Ax ~ 
a/ kBTAt/(a(±). We were surprised to find instead that there was no measurable deviation 
of P(r) from the theoretical predictions, even at extremely small values of r, as shown in 
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Figure El Our understanding of this result is based upon an interpretation of the algorithm 
as a form of Monte Carlo sampling. 

B. Monte Carlo 

Monte Carlo simulation is a method of efficiently sampling a desired equilibrium proba- 
bility distribution P eq (X) for some set of coordinates X. A Monte-Carlo simulation is based 
on a Markov process in which the conditional probability T(Xi — > Xf) of a transition from 
a state Xi to a state Xf is related to the probability of the reverse transition by a detailed 
balance condition [3] 

P eg (X i )T(X i -> X f ) = P eq (X f )T(X f -> X«) • (25) 

If the Monte-Carlo jump process involves both the generation of a trial move, and a decision 
regarding whether to accept or reject the move, then 

T{Xi -> Xf) = G(Xi -> -> X/) (26) 

where G(Xi — > X/) is the probability of generating a particular trial move and — > Xf) 
is an acceptance probability. 

1. Uncrossable Threads 

We first consider the algorithm used to simulate infinitely thin uncrossable polymers. The 
only potential energy in the problem is the intramolecular bending energy Uq of all of the 
chains. In our algorithm, trial moves are generated by a Brownian dynamics algorithm that 
has been designed so as to generate an equilibrium distribution for this potential. We thus 
consider a class of algorithms in which a BD algorithm is used to generate trial moves, and 
in which the BD algorithm is assumed (for the moment) to approximately satisfy detailed 
balance in the limit At — > 0. That is, we assume for the moment that 

PoiXjGiXi - Xf) = P (Xf)G(X f - Xi) (27) 

in the limit At — ► in which the BD algorithm is designed to yield the equilibrium dis- 
tribution P (X). Here, Pq(X) oc e~ u °( x ^ is the equilibrium distribution for a solution of 
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non-interacting wormlike chains, which is also the correct equilibrium distribution for a 
solution of infinitely thin but uncrossable chains with no interaction potential. 

The prohibition on chain crossing is enforced by simply rejecting all moves that cause 
one chain to cut through another (A = 0) and accepting all others (A = 1). If a transition 
Xi — > Xf involving the motion of the beads of one chain causes that chain to cut through 
another, then the hypothetical reverse move Xf —>■ Xi, which would involve moving the 
same chain from its final to initial state, would also cause a chain crossing, and so would 
also be rejected. If the matrix G satisfies detailed balance, but certain types of moves 
are prohibited, then the overall transition matrix T = GA will thus also satisfy detailed 
balance, as long as the reverse of any prohibited move is also prohibited, as is the case here. 
An algorithm in which a BD algorithm is supplemented by such a rejection scheme will 
thus satisfy detailed balance to the same extent as the underlying BD algorithm. If the BD 
simulation can be shown to satisfy (or nearly satisfy) detailed balance, this would explain 
the otherwise surprising accuracy of our results for P(r). 



2. Intermolecular Interactions 



Our algorithm for simulating a system of chains with a short range repulsive intermolec- 
ular repulsion is an extension of the above idea. We express the total potential energy of 
the system as a sum 

U = U + U int (28) 

where Uq is a relatively "soft" intramolecular bending energy and Ui n t is a "hard" inter- 
molecular interaction energy. The equilibrium distribution for the interacting system is of 
the form 

P eq (X) = P (X)e- u ^ )/kBT (29) 

where Pq(X) is the equilibrium distribution for the non-interacting gas, which is governed 
by the bending energy. If moves are generated with a BD algorithm that satisfies detailed 
balance condition (|27j) for non-interacting polymers, then the corresponding condition for 
the interacting solution can be satisfied by requiring the acceptance probability A to satisfy 

A(Xj -» X f ) _ e _Ac/ int /fe s r (30) 
A(X f ^XA 
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where AUi nt = Ui nt (Xf) — Ui nt (Xi) is the change in interaction energy. We accomplish 
this by supplementing our prohibition on chain crossing by the usual Metropolis scheme of 
taking A(X { -»• X f ) = 1 for AU mt < and A(X { -> X f ) = e -*u int /k B T for A[/ .^ > Q for 
moves that do not cause chains to cross. The same reasoning underlies our slithering snake 
algorithm for equilibrating such systems, in which the trial moves are slithering snake moves 
that are designed to satisfy detailed balance for non-interacting chains. 

Our use of a Brownian dynamics algorithm to generate trial moves and a Metropolis 
acceptance scheme based on the interaction energy is useful in situations in which we wish 
to resolve the dynamics of chain bending, and At can be made small enough to do so, but 
in which we do not need to resolve the dynamics of intra-chain collisions, and in which the 
range of repulsion between chains is so small that dynamically resolving collisions would 
require us to use a much smaller time step. Analogous hybrid algorithms are potentially 
useful in any problem in which one wishes to simulate diffusion in a comparatively soft 
potential in the presence of an additional steeply repulsive potential. 

3. Do BD simulations satisfy detailed balance? 

In appendix [BJ we try to answer the question of whether a valid Brownian dynamics 
algorithm satisfies detailed balance by considering a simplified model of one dimensional 
(ID) diffusion under the influence of a soft potential Uq{X). There we consider both an 
algorithm in which the random force is chosen from a uniform distribution and an algorithm 
in which the random force is chosen from a Gaussian distribution. We find that, for this 
model, both algorithms satisfy detailed balance for simple diffusion in a constant potential, 
with no drift. In problems with nonzero deterministic drift forces, it is straightforward 
to show that the BD algorithm with normally distributed random forces satisfies detailed 
balance for any linear potential, and should thus satisfy detailed balance in the limit of 
At — > for any locally linear potential. It is equally easy to show that an algorithm with 
a uniformly distributed random force in addition to a deterministic force does not exactly 
satisfy detailed balance. It is argued, however, that the small violation of detailed balance 
found with uniformly distributed noise should have a negligible effect when the step size 
is taken to be small enough to resolve variations in the soft potential. We thus conclude 
that, with either type of BD algorithm, the effect of any violation of detailed balance should 
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vanish in the limit At — > of interest. 

Our conclusions about the simple ID diffusion problem are consistent with our results 
for simulations of entangled polymers, in which we find that the combination of single- 
chain BD simulation with an appropriate rejection criterion leads to results that agree to 
within stochastic error with theoretical predictions for the radial distribution and other 
equilibrium properties. In most of our simulations, we have used an unprojected random 
force r\' whose Cartesian components are chosen from a uniform distribution. An extensive 
series of simulations was carried out with this choice before we understood its potential 
theoretical significance. In order to check whether our use of uniform distribution affected 
our results, we also ran a smaller number of simulations with a Gaussian distribution for 
r) 1 . Figures |H1 and El compare results for P(r) and (Ad 2 (t)) obtained from simulations of 
infinitely thin threads using these two different noise distributions. For the value of At used 
in our simulations, the results agree to within statistical error. Because the mathematical 
argument required to show that detailed balance is obtained in the limit At — > is relatively 
straightforward only for BD algorithms with Gaussian-distributed random forces, however, 
we recommend that any future applications of a hybrid BD-MC algorithm analogous to 
that used here use Gaussian distributed random forces, if only to simplify any subsequent 
discussions of the underlying theory. 



VI. CONCLUSIONS 



We have presented a hybrid algorithm for simulating highly entangled solutions of very 
thin uncrossable wormlike filaments. The algorithm uses a single-chain Brownian dynamics 
(BD) algorithm for generating trial moves. This is supplemented with a Monte Carlo-like 
scheme where any trial move that causes two chains to cross or that leads to excluded volume 
overlaps are rejected. Efficient methods were developed for checking for chain crossings and 
excluded volume violations. Because the algorithm is based on a Brownian dynamics time- 
stepping algorithm, it can be used to simulate dynamical, as well as equilibrium properties. 
The algorithm is expected to be much more efficient than Brownian dynamics simulation 
of a bead-spring model of nearly tangent beads in the limit of highly entangled solutions of 
extremely thin chains. Comparison with the efficiency of a molecular dynamics simulation of 
a bead spring model is complicated by the difference between the diffusive dynamics of our 
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FIG. 8: Effect of noise distribution on P(r). cL 3 = 1000, iV = 20, L p = L, (\\ = C±, At/r"^ = 
4.5 x 10" 8 . 

algorithm and the ballistic motion obtained in molecular dynamics of a semi-dilute solution 
in the absence of explicit solvent, which would lead to efficient sampling of configuration 
space but a loss of dynamical realism. We have used the algorithm discussed here to study 
dynamics and stress relaxation in highly entangled solutions ^(|, which will be discussed 
elsewhere. 

APPENDIX A: RADIAL DISTRIBUTION FUNCTION FOR RODS 

Consider a test rod oriented along the z-axis. Let r be the closest approach vector and 
n be the unit vector normal to the plane defined by the test rod and the closest approach 
vector. Choose an infinitesimal area on this plane, dA = dr x dz. Since dr is perpendicular 
to dz and n is perpendicular to this area unit, one can write dA = drdzn. Using geometry, 
one can show that the number of rods, with orientations lying between u and u + du, that 
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FIG. 9: Effect of noise distribution on (Ad 2 (t)}. cL 3 = 1000, N = 20, L p = L, C]| = (±, 
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pierce this infinitesimal area is given by 



dN u = p dA ■ u (2) 



Arc ' 



(Al) 



where p, as defined earlier, represents the contour length density. d 2 u represents the solid 



angle subtended by u and u + du and is given by d 2 u = sin 



, where 9 represents 



the angle between u and the z-axis and represents the azimuthal angle. The factor 4ir 
in the denominator is obtained by integrating with respect to d 2 u over the surface of a 
sphere. The factor of 2 in the numerator arises because of the need to count for rods with 
both orientations, u and — u. Substituting for d 2 u and dA from above while noting that 
n ■ u = cos(90° — 6) = sin8 and integrating equation (jAljl over dd, d<j), we obtain the radial 
distribution function given in Eq. (|2"Tf . 
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APPENDIX B: SIMULATION OF ONE-DIMENSIONAL DIFFUSION 



In this appendix, we consider whether a very simple Brownian dynamics algorithm satis- 
fies detailed balance when applied to one dimensional (ID) diffusion. We consider the use of 
a hybrid BD-MC algorithm analogous to that discussed in the body of the paper to describe 
the diffusion of a particle with a position X(t) which is subjected to both a slowly-varying 
(or "soft") potential Uq(X) and an arbitrary steep repulsive (or "hard") potential U\(X) 
that represents the walls of a confining ID box. To represent a box of length L with hard 
walls, we might take U 1 (X) to vanish for — L/2 < X < L/2, and to diverge outside this do- 
main. We consider an algorithm in which trial moves are generated by a Brownian dynamics 
algorithm that is obtained by discretizing the Langevin equation 

(dX= [F (X)+rj(t)]dt (Bl) 

where ( is a friction coefficient, X is the particle position, rj(t) is a random force, and 
F (X) = —dU(X)/dX is a force arising from the soft potential. Each time step, we generate 
a trial step Xi Xf of magnitude 

(AX = [F (X t )+r ] ]At (B2) 

where AX = Xf — JQ, and where 77 is chosen from a distribution P(rj) with moments 
(77) = and {if) = 2ksT(/At. Here, we consider whether the resulting jump probability 
satisfies the detailed balance condition P (X i )G(X i — > Xf) = P (Xf)G(Xf — > JQ), where 
Po(X) oc e - u o( x )/ k sT _ jf ; j n a given step, we generate a trial move Xi — > Xf using a force 
we must consider a hypothetical reverse move in which the random force r]~ must satisfy 



r ] + + F (X l ) = -[r ] - + F (Xf)] (B3) 

in order to generate a displacement —AX. 

If the probability distribution P(i]) for the random force 77 is a Gaussian, 

P(rj) oc exp (-r) 2 /(4k B T()) , (B4) 

it is straightforward to show that 

GjXj -> Xf) = e -F[ v +F]At/(Ck B T) (g 5 ) 

G(X f X^ 
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where F = [F (Xj) + F (Xf)]/2. To lowest order in AX, we may approximate the change 
in energy AU = U(X F ) - U(Xi) by AU ~ F5X and approximate AX ~ [77 + F]At to show 
that 

G(Xj -» X/) _ AC/o/fc s T / B6 n 

The algorithm thus approximately satisfies the detailed balance condition in the limit At — > 
of interest, and can be shown to exactly satisfy detailed balance in the special case of a 
constant force F (X). 

We next consider a uniformly distributed random force, for which P(rj) = 1/(2E) over a 
range — E < rj < E, where E = ^Q^bT / At. This algorithm clearly does not exactly satisfy 
detailed balance, since a step for which the random force r/ + that generates the forward trial 
move is near one of the boundaries ±E can sometimes be reversed only by a random force r/~ 
that lies outside the allowed range [-E, E\. An algorithm in which the reverse of an allowed 
step is sometimes not allowed clearly cannot exactly satisfy detailed balance. In the limit 
of small force F , however, this violation becomes small, insofar as it only affects a small 
fraction of all jumps for which r] falls within a narrow range of values near its maximum or 
minimum allowed value of ±E. In the special case of a vanishing force, or constant potential 
Uo(X), detailed balance is recovered, for any distribution with the symmetry P{rj) = P(—r)), 
including a uniform distribution. 

In order to assess the importance of this violation of detailed balance for an algorithm 
with uniformly distributed noise and a spatially varying potential Uq, we must consider 
the implications of our original assumption that Uq(X) is a "soft" potential and U±(x) is 
"hard". The effect of a deterministic force F becomes comparable to or greater than that 
of diffusion only when applied over a characteristic distance of order ksT/F or greater, for 
which the corresponding potential energy exceeds ksT, or over a corresponding time scale 
of order ksTC/F 2 , for which the displacement Ft/( caused by drift exceeds the root-mean- 
squared displacement \f~Dt caused by diffusion. If the range of the hard potential Z7i is 
much less than the characteristic length scale ksT/F, then we expect the existence of a 
nonzero force F to have negligible effect upon the structure of the probability distribution 
within the narrow boundary layer in which the hard repulsive potential is important, other 
than to change an overall prefactor that is sensitive to the nature of the far field. Within 
this boundary layer, the form of the probability distribution is instead determined by a 
balance between the repulsive potential U\ and diffusion. Since the form of the equilibrium 
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distribution produced by this balance can be correctly simulated by an algorithm with no 
drift, which does satisfy detailed balance, we expect the addition of a small drift term (and 
a correspondingly small violation of detailed balance) to have little effect upon the form of 
the probability distribution near the boundary, as long as the range of the "hard" potential 
U\ is much less than the characteristic length scale ksT/F of the "soft" potential. 

For completeness we may also consider the case in which the range of the repulsive 
potential U 1 is actually comparable to that the characteristic length scale of the soft potential 
Uq. Since we rely on a BD algorithm to resolve the effect of potential U , we assume that 
we will choose a step size At such that the corresponding spatial step AX ~ V DAt is much 
less than any length characteristic of potential Uq. If the characteristic length scales of Uq 
and U\ are comparable, however, then we will automatically have chosen the step size small 
enough so that potential Ui changes by much less than ksT per step. It can be shown in this 
case, by considering the effect of the acceptance criterion in the limit of a soft potential U±, 
that the overall transition matrix T = GA for the resulting algorithm satisfies the conditions 
on the first and second moments necessary to obtain a valid Brownian dynamics algorithm 
for the total potential Uq + U\. In this case, the hybrid algorithm can thus be justified as a 
valid form of Brownian dynamics simulation, and the violation of detailed balance becomes 
irrelevant. 



30 



[1] K. Kremer and G. Grest, J. Chem. Phys. 92, 5057 (1990). 

[2] R. Everaers, F. Julicher, A. Ajdari, and A. C. Maggs, Phys. Rev. Lett. 82, 3717 (1999). 

[3] M. Pasquali, V. Shankar, and D. C. Morse, Phys. Rev. E 64, 020802 (2001). 

[4] P. Dimitrakopoulos, J. F. Brady and Z. G. Wang, Phys. Rev. E. 64, 050803(R) (2001). 

[5] V. Shankar, M. Pasquali, and D. C. Morse, J. Rheol. 46, 1111 (2002). 

[6] M. Pasquali and D. C. Morse, J. Chem. Phys. 116, 1834 (2002). 

[7] A. Montesi, M. Pasquali, and D. C. Morse, J. Chem. Phys. 122, 084903 (2005). 

[8] E. J. Hinch, J. Fluid. Mech. 271, 219 (1994). 

[9] M. Fixman, J. Chem. Phys. 69, 1527 (1978). 

[10] D. C. Morse, Advances in Chemical Physics 128, 65 (2004). 

[11] M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Clarendon Press, 1996). 

[12] G. Lamm and K. Schulten, J. Chem. Phys. 78, 2713 (1983). 

[13] H.-C. Ottinger, J. Chem. Phys. 91, 6455 (1989). 

[14] E. Peters, Phys. Rev. E 66, 056701 (2002). 

[15] K. Binder, Monte Carlo and Molecular Dynamics Simulations in Polymer Science (Oxford 

University Press Inc., 1995). 

[16] S. Ramanathan, Ph.D. thesis, University of Minnesota (2006). 



31 



